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ABSTRACT 

In numerical models of thin astrophysical disks that use an Eulerian scheme, gas orbits supersonically 
through a fixed grid. As a result the time step is sharply limited by the Courant condition. Also, 
because the mean flow speed with respect to the grid varies with position, the truncation error varies 
systematically with position. For hydrodynamic (unmagnetized) disks an algorithm called FARGO 
has been developed that advects the gas along its mean orbit using a separate interpolation substep. 
This relaxes the constraint imposed by the Courant condition, which now depends only on the peculiar 
velocity of the gas, and results in a truncation error that is more nearly independent of position. This 
paper describes a FARGO-like algorithm suitable for evolving magnetized disks. Our method is second 
order accurate on a smooth flow and preserves V • B = to machine precision. The main restriction is 
that B must be discretized on a staggered mesh. We give a detailed description of an implementation of 
the code and demonstrate that it produces the expected results on linear and nonlinear problems. We 
also point out how the scheme might be generalized to make the integration of other supersonic/super- 
fast flows more efficient. Although our scheme reduces the variation of truncation error with position, 
it does not eliminate it. We show that the residual position dependence leads to characteristic radial 
variations in the density over long integrations. 
Subject headings: numerical methods, magnetohydrodynamics 



1. INTRODUCTION 

Numerical experiments have played a key role in ad- 
vancing our understanding of accretion disk dynamics. 
Future attacks on the main unsolved problems of disk 
theory, such as the evolution of large-scale magnetic fields 
in disks, will also likely benefit from numerical experi- 
ments. But numerical work is always limited by current 
hardware and algorithms. Here we describe a new al- 
gorithm for evolving the magnetohydrodynamic (MHD) 
equations that is designed to speed up, and improve the 
quality of, future disk experiments. 

The majority of numerical hydrodynamic studies of 
disks use an Eulerian approach: the fluid equations are 
discretized in a fixed frame and the code "pushes" the 
fluid through the grid. In an accretion disk the fluid 
velocity can be written 

V = Vorb + At) (1) 

where Vorb is the circular orbit velocity, and Av rep- 
resents departures from a circular orbit caused by, e.g., 
turbulence. If the equations are discretized in a nonro- 
tating frame then typically an Eulerian scheme will have 
to push fluid through the grid with a speed that varies 
systematically with radius r. As a result the truncation 
error will vary with position, possibly yielding misleading 
results. 

Another disadvantage of a straightforward Eulerian ap- 
proach has to do with the Courant condition on the time 
step At: 

At < (2) 
VAL ^ ' 

where C ~ 1 is the Courant number, V is the fastest wave 
speed in the problem, and AL is the grid scale. If one 
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is interested in a cold accretion disk with \Av\ ^ Cs ^ 
\vorb\ [cs = sound speed), then V will be dominated by 
the orbital motion, i.e. V ~ rfl. This implies small time 
steps: the code is only allowed to push the fluid along its 
orbit by a fraction of a zone per time step. Most of the 
computational time will be spent on orbital advection. 

All this runs contrary to one's sense that somehow the 
peculiar motion Av should control the time step and 
the truncation error. After all, in a frame moving on a 
circular orbit one expects that the motion of the fluid is 
either subsonic or near-sonic. 

Three strategies have been employed to get around 
these two issues. First, one can work in a rotating frame. 
This works well only within a few scale heights H = Cs/^ 
of the corotation radius. Second, one can employ a 
Lagrangian scheme, which follows individual fluid ele- 
ments. In astrophysical applications this usually means 
smoothed particle hydrodynamics (SPH). SPH has in- 
trinsic noise that makes it unsuitable for many sensitive 
disk dynamics problems. It is also difficult, in our expe- 
rience, to incorporate magnetic fields into SPH. Third, 
one can adopt a hybrid, quasi-Lagrangian scheme that 
treats the orbital a dvec tion separately. This is the ap- 
proach adva nced by|M assct (200 3) in his FARGO code , 
and later by Gammic (2001), Joh nson fc Gammid ()2003D 
and Johnson & Gammic (2005). Our contribution here 
is to extend this method to MHD. 

There are two other possible strategies that we are 
aware of for addressing these problems in the context of 
disks. First, one can do an orbitally-centered domain 
decom position in a parallelized code (Gaunt fc K o r"pi| 
1200 If ). Each processor works on a small portion of the 
grid (< H) without orbital advection, and then orbital 
advection is used in the boundary conditions to link the 
small portions together. Second, one could employ a 
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fully Lagrangian orbital advection by defining the grid 
in shearing coordinates, coupled to a remap (similar to 
what is done here) once per shear time {qft)~^ (Narayan, 
private communication) at the risk of introducing a new, 
numerical timescale into the problem. This is what is 
done in so me spectral schemes for i ncompressible shear 
flows (e.g.. lUmurhan fc Regevll2004D . 

The main idea in the approach we take here is to 
operator-split the update of the fluid variables. The evo- 
lution equation for each dependent fluid variable F is 



dF 

^^-(.•V)F + A 



(3) 



where the first term on the rhs is advection (its form is 
determined by Galilean invariance) , and the second term 
C contains everything else. The advection operator can, 
in turn, be split again: 

^ = -{vorb-V)F~{AvV)F + C, (4) 
ot 

The orbital advection operator simply pushes the fluid 
elements along its orbit. Remarkably, this can be done 
using an interpolation formula that is not constrained by 
the Courant condition! One simply needs to know 



x{t) 



dt'Vorb{t') 



(5) 



in advance, so that the fluid variable F{x[t -\- At]) can 
be interpolated from known values near F{x[t]). The 
method can be made formally second-order accurate us- 
ing Strang splitting. Notice that this idea can be applied 
to any flow with known orbits, not just circular, Keple- 
rian orbits for disks. 

Implementing a stable, accurate orbital advection op- 
erator involves a surprising amount of bookkeeping, par- 
ticularly when one must maintain a divergence-free mag- 
netic field. In this paper we describe an implementation 
for a particular context, albeit one of considerable in- 
terest: the "local model" for astrophysical disks. The 
plan of the paper is as follows. In §2 we write down 
the basic equations describing our model, and show how 
the advection can be split into pieces corresponding to 
the orbital and peculiar velocities. In §3 we summarize 
our algorithm, deferring its rather tedious derivation to 
Appendix |^ In §4 we describe tests of the method. §5 
describes a sample nonlinear application. We have incor- 
porated our algorithm into a ZEUS-like code for perform- 
ing calculations. §6 summarizes the results and identifies 
the key formulae for implementing our scheme. 

2. BASIC EQUATIONS 

The "local model" for astrophysical disks is obtained 
by expanding the equations of motion around a circular- 
orbiting coordinate origin at cylindrical coordinates 
(r, (p, z) = (ro, f2ot-|-(/)o, 0), assuming that the peculiar ve- 
locities are comparable to the sound speed and that the 
sound speed is small compared to the orbital speed. The 
local Cartesian coordinates are obtained from cylindrical 
coordinates via (cc, y, z) = (r — Tq, ro[4> — ^ot — (j)o], z). 

In this context the equations of isothermal ideal MHD 
consist of seven evolution equations, given by 



dp 
dt 



V • (pv) = 0, 



(6) 



dv 

'dt' 



,Vp {B-V)B 



p Slip 



Anp 



dB 

'dt 



V X {v X B) 



(7) 

0, (8) 



plus the divergence-free constraint on the magnetic field: 

V ■B = 0. (9) 

The final two terms in equation ([7]) represent the Coriolis 
and tidal forces in the local frame. The orbital velocity 
is 

Vorb = -Q^xy, (10) 

where 

is the shear parameter. One can readily verify that this 
velocity, along with a constant density and zero magnetic 
field, is a steady-state solution to equation ([7]). 

Integrating equation ([9]) over a control volume and ex- 
pressing the volume integral as a surface integral via 
Gauss's Law gives an alternative representation of the 
divergence- free constraint: 



$ = / B • n da = 0, 

JK 



(12) 



where A is the surface bounding the volume, da is an 
area element in that surface and n is a unit vector nor- 
mal to the surface. Satisfying expression (fT2|) throughout 
the evolution of equations ©-([S]) is one of the main chal- 
lenges in numerical MHD. 

The evolution equations ([6])- ([8]) can be recast using 
V = Vorb + At): 



dp 
dt 



Vp + V ■ (pAv) ^ 0, 



(13) 



dt 

VB^ B ■ VB 



Btt/O Anp 



P 



2nx Av- qn (Av)^ y = 0, (14) 



Vorb ■ -V x{Avx B)+ qilB,^ y = 0, (15) 



dB 

'dt 

There are three differences between equations (|T3| - (|15p 
and the original equations (O-®: 1) each equation has 
an additional transport term due to the orbital (mean 
shear) velocity, Vorb ■ V, 2) the tidal term in equation 
([7]) has been replaced by —qil (Av)^ y in equation (|14p 
and 3) there is an additional term q^B^y in equation 
(fT5|) . The latter two terms reflect the conversion of radial 
velocity and magnetic field components into azimuthal 
components by the shear. The last term in equation (jl4p 
can simply be treated as an additional term in the finite- 
difference algorithm, whereas the last term in equation 
(|15p must be treated differently, using the algorithm we 
outline in this paper, in order to preserve the divergence- 
free constraint. 

The local model is usually simula ted using the "shear - 
ing box" boundary conditions (e.g. iHawlev et al.l ll995'). 
These boundary conditions isolate a rectangular region 
in the disk. The azimuthal [y] boundary conditions are 
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periodic; the radial (x) boundary conditions are "nearly 
periodic" , i.e. they connect the radial boundaries in a 
time-dependent way that enforces the mean shear flow; 
and one is free to choose the vertical boundary conditions 
for physical and numerical convenience. 

3. ALGORITHM 

The orbital advection substep consists of evaluating 
the fluid variable F at time t + At using interpolation: 

F{x, y,z,t + At) ^F{x,y + qilxAt, z, t) . (16) 

Methods for stable interpolation of the independent vari- 
ables are well known. An implementation for hydro- 
dynamic a l vari ables is described for the local model by 
iGammid (l200ll) . 

The interpolation can be thought of as shifting a single 
column of zones (at constant x and z) by what is gen- 
erally a noninteger number of zones. The shift can then 
be decomposed into an integer number of zones and a 
fractional shift of up to half a zone. The integral shift 
can be done trivially, while the fractional zone shift is 
best done using the same transport algorithm as the rest 
of the code. 

The induction equation must be treated differently be- 
cause of the V • B = constraint. In our code the mag- 
netic field is discretized on a staggered mesh, and mag- 
netic field variables represent fluxes through zone faces. 
The effect of orbital advection on the zone faces in the 
X — y plane is illustrated in Figure [TJ The dashed lines 
show the positions of the "old" zone faces after they have 
been sheared through a time At. The solid lines show 
the "new" zone faces onto which the fluxes from the old 
zone faces must be interpolated. Flux freezing requires 
that the fluxes through the old zone faces be preserved 
by the orbital advection; our algorithm simply interpo- 
lates the fluxes in these sheared zones onto the new zones 
in a way that preserves V • i? = 0. 

3.1. Definitions 

The shear has two effects on each zone of the old grid: 
1) a linear distortion of the zone in the azimuthal di- 
rection, and 2) an azimuthal advection of the zone that 
depends upon the radial position of the zone in the old 
grid. We quantify these two effects with the following 
definitions: 

q^AxAt 

s = 1 (17 

Ay ^ ' 

is the relative shift (in dimensionless zone units) of a fluid 
element across a single zone in one time step {—sAy/ Ax 
is the slope of the diagonal fines in Figure [J), and 

VorbAt -qilxAt 

is the amount (in dimensionless zone units) that a fluid 
element is advected by the shear in one time step. In 
general, S is composed of a non-integral number of zones, 
which we divide into an integral part NINT(5') and a 
fractional part 

f = S- NINT(S'). (19) 

Here NINT(S') is the value of S rounded to the nearest 
integer, so that / can take on both positive and negative 
values. 



We denote old zones by the superscript n and new 
zones by the superscript n + 1. The indices i,j,k corre- 
spond to the x,y,z directions. The azimuthal index for 
an old zone goes from j to 

J = j- NINT(S') (20) 

after each time step At. 

3.2. Interpolation Formulae 

We can obtain divergence-free interpolation formulae 
by considering a control volume (which we will call a 
subvolume., because it is smaller than a zone) bounded 
by portions of zone faces from both the old grid and the 
new grid (which we will call subfaces, because they are 
in general smaller than a full zone face), and requiring 
that the sum of the fluxes into or out of that volume is 
zero. 

There are three distinct cases that occur when mapping 
the old, sheared grid onto the new grid, depending on 
whether / is positive or negative and whether or not 
the azimuthal face of a sheared grid zone intersects the 
azimuthal face of a new grid zone. The three cases are 
illustrated in FiguresHHH and correspond to |/|/s < 1/2 
(Case 1), f/s > 1/2 (Case 2) and f/s < -1/2 (Case 3). 
The value of //s depends, in turn, on the x coordinate 
of the zone and the time step (see Figure [1]) . 

Deducing the fluxes through the faces of the new 
zones is a matter of frankly tedious bookkeeping that 
is described in detail in Appendix |X] but summarized 
here. First, in each case write the constraint that the 
sum of the fluxes in and out of each subvolume vanish 
(V ■ B — 0). Next, solve for the unknown fluxes through 
the subfaces of the new grid in terms of the fluxes through 
the subfaces of the old grid. The latter can be deduced 
given a model for the variation of the field strength over 
each zone face in the old grid; we use a linear model with 
van Leer slopes, consistent with the rest of our ZEUS- 
like code. Finally, sum the fluxes through the subfaces 
to obtain the flux through each face of the new grid. 

The final formulae are given below, where the w 
coefficients appearing in these expressions are defined in 
Table Al and the dq's are van Leer slopes: 

Radial field update. Case 1: 

b^7ji' = Kjfe + ^io(&<j-ife - bx-j,) 

+ wii{hx_dqy'lj_^k ^ bx.dqylj^), (21) 

Radial field update. Case 2: 

K^k = KVfc + i«2o(Kj-ife - bxljk) 

+ W2i{bx.dqy2j_ik - bxAqylj^^), (22) 

Radial field update. Case 3: 

b^'itk ^ bxljk + w^3o(&a;"j+ife - hxl^^,) 

+ W31 {hxAqylj^^^ - bx.dqy^jf^). (23) 
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Azimuthal field update, Case 1, n even: 



by, 



n+l 
ijk 



1 



Ay r 



{bz^jk+i bz^j^) 



Ay 



+ Fi4(K"/-i/c+i - bzZ-ik) + wnibz^jk - bz^jf,^^) 

+ wi^{hz_dqxlj_^k+i - bz-dqx"j_j^f,) + wi(i{hz_dqylj_^k+i 
-bz-dqylj_ik) + Wi8{bz.dqXij,, - bz_dqxljk+i) 
+ wig{bz.dqyijk - bz.dqy1jk+i)] ■ (24) 



Azimuthal field update, Case 1, n odd: 



by'-'t' = - 



by?jk + by: 



wubx. 



i+ljk 



Wiibx-dqy, 



Ay r 

— [-wioKz-ife 
Wi3bx_dqy^^i 



ij-lk 



Ay 
Az 



Wi4{bz, 



ij-lk+l 



K"/-ife) + wi5(bz_dqxlj_j_^^^ 



- bz_dqx'^j_^k) + wi6(^'z_dqy,"j_ifc+i - bz_dqy'^j_^^.) 
+ wi7(K"jfc - Kjfc+i) + wisibz.dqx'ljk - bz.dqx'lji,^-^) 
+ wigibz.dqyljk - bz.dqy^jk+i)] . (25) 



Azimuthal field update. Case 2: 



by 



n+l 
ijk 



1 

2 
Ay 
Az 



byuk + &2/z"j-ifc + ^ (KV-ifc - 6a;r+ij-ifc) 



(b^ij-ik bz^j ^ii^j^-i) 



Ay 



~W2obXi 



ij-lk 



+ M;22&a;f+ij_ife - W2ibx_dqylj_^^ + W23bx.dqy^_^^j_^^] 
Ay 

+ [w24{bz^j-ik+i - bzlj_^^) + u;25(&^-rfg</_ife+i 
bz-dqxlj_-^f.) + W26{bz.dqy2j _ik+i - bz.dqylj_-^k)\ ■ (26) 



Azimuthal field update. Case 3: 



byljk + &^/^"J+lfc + ^ {bx'l+ijk - bxfjk) 



Ay r 

— [w30teZ/fc - W32bx^_^_i 



Jk 



+ wi^{bzjlqxlj_^k - bz.dqxljf.) + wifi{bzjlqy'^j_-^f, 
- bzjlqyljt,) + Wls.{bz.dqxlJ+^^. - bzjiqx^j^) 

+ wigibz.dqy^j^^^ - bz.dqy^j^), (28) 

Vertical field update, Case 2: 



^^i'jk^ — bZuk 



W24(64j-ifc - bzijk) 



+ W25{bz_dqx'lj_^k - bz-dqx2jk) 
+ W2Q{bz.dqylj_-^k - bz.dqylji,), 



(29) 



%jk 



Vertical field update, Case 3: 

= bzlji, + u.34(K" +ife - bzljk) 
'W3^{bz_dqxlj^^k - bz_dqxijf,) 
-W3Q{bz.dqylj^^^ - bz.dqylj^). 

4. TESTS 



(30) 



Ay 

wsibx.dqy'^jk - u;33&a;-rfgyIVi,/fc] ^ [w34{bz^jk 
-bzl,f,^^)w3<i{bz_dqxljk - 62_dga;"jfc+i) 
+ {bz.dqyljk - bz.dqyljk+i )] . (27) 



Vertical field update, Case 1: 
K-fc' = ^^^"^fc + "'i4(K'j-ifc - + w'i7(K"/+ifc - bz^jk) 



We have tested our algorithm on both linear and 
nonlinear problems. Linear perturbations in the local 
model are decomposed most naturally in terms of shear- 
ing waves, or shwaves, which appear spatially as plane 
waves in a frame comoving with the shear. The radial 
wavenumber of a shwave increases linearly with time and 
its amplitude does not in general have an exponential 
time dependence (as does a normal mo de). D etails on 
shwaves in isothermal MHD are given in [Johns on (200^ 
and summarized in Appendix |B1 We have calculated the 
evolution of both compressive and incompressive shwaves 
as a function of numerical resolution, and the results are 
shown in Figures [MTTl 

We employ a grid of physical size x 
and numerical resolution x Ny x N^. The equi- 
librium state about which we perturb has a constant 
density po and spatially constant magnetic field Bq, 
plus the background shear flow. On this background 
state we impose a plane wave perturbation with ini- 
tial amplitude (5p,5v,5B) and initial wavevector k = 
2'n:{mx/ Lx,my / Ly,mz/ Lz), where mx/iriy < corre- 
sponds to a shwave that initially leads the mean shear. 
The perturbations are expressed in units with ~ Cg ~ 
1. 

Our first linear test is a simple advection of the mag- 
netic field components with zero velocity perturbation, 
for a shwave that swings from leading to trailing. Differ- 
ent operators in an operator split scheme do not neces- 
sarily converge at the same rate; the overall convergence 
rate depends upon the combined convergence properties 
of each operation. This test is therefore important for 
isolating the convergence properties of our algorithm. 
For this test, we employ equal box dimensions L = AH 
and equal numerical resolutions N . The other param- 
eters for this run are rrix = —1, niy = rn^ = 1, and 
i?o = 0. The initial perturbation in 5p — 5v — Q and 
5B = 10^6(2, l,l)cos(fc • I), where I = {x,y,z). The 
amplitude of these shwaves is constant with time. Fig- 
ure [5] shows the evolution of the vertical field component 
at iV = 8, 16, 32 and 64. 

The convergence properties of our algorithm for this 
test are shown in Figure [U which is a plot of the LI 
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norm of the error in each magnetic field component as 
a function of numerical resolution N. Also shown on 
this plot are the convergence properties of a run with 
orbital advection turned off, for comparison purposes. 
The algorithm converges at second order, as expected. 

To demonstrate the improved accuracy obtained by us- 
ing orbital advection, we have run the same test at vari- 
ous box sizes. Figure [7] shows the LI norm of the error in 
the azimuthal field component in runs with L — H and 
L = lOH, both with and without orbital advection. The 
errors are comparable in the run with L = H , but in the 
run with L — lOH the errors with orbital advection are 
smaller by a factor of ^ 4. For our second-order algo- 
rithm, this corresponds to a gain in effective resolution 
(at fixed error) of ~ 2 (for L = lOH). Orbital advection 
is more efficient in addition to being more accurate, par- 
ticularly when the box size is large compared to H . For 
example, at iV = 64, the ratio of zone cycles with orbital 
advection on and off is ~ 0.8 in runs with L = H; this 
ratio decreases to ~ 0.2 in nms with L = lOH. 

Figure [5] shows the evolution of the radial field pertur- 
bation for an incompressive shwave that grows nearly 
exponentially as it swings from leading to trailing. 
The parameters for this run are — Ly — lOH, 
Lz ^ H, rrix ^ -2 , rUy = ^ I, Ny = = 
iV^/2, and Bq = ^Jl^/lQ[n/k^)z} The initial per- 
turbation is 8p = 8.95250 x 10"i''cos(fc • I - 7r/4), 
8v = 10-8(8.16589, 8.70641, 0.762537) cos( A; • I + tt/4), 
and SB = 10-^(-1.08076, 1.04172, -0.320324) cos(fc • I - 
7r/4). 

As discussed by iJohnson fc Gammii (|2005l ). aliasing 
of incompressive shwaves can artificially convert trail- 
ing shwaves into leading shwaves. Figure [9] shows the 
long term evolution of the previous run, demonstrating 
that aliasing can result in artificial growth in the linear 
regime. We do not consider this to be a serious problem 
for a nonlinear calculation, however, such as the devel- 
opment of turbulence due to the MRI. The growth rate 
due to aliasing cannot exceed the MRI growth rate in the 
linear regime, and the evolution in the nonlinear regime 
is dominated by small scale fiuctuations that interact on 
a time scale much shorter than the shear time scale. In 
addition, the strong aliasing seen in Figure [9] depends 
upon the very small amount of diffusion present in this 
test due to the lack of any motion with respect to the 
grid. To introduce numerical diffusion, we perform the 
same test with an additional bulk epicyclic motion of 
the grid superimposed (amplitude O.lcg). As shown in 
Figure 1101 a small amount of diffusion can significantly 
reduce the effects of aliasing. 

Figure [Tl] shows the evolution of the azimuthal field 
perturbation for a compressive shwave. The parameters 
for this run are L ~ 0.5H, nix = —2, my = = 
1, Ny ^ ^ Nx/2, and Bo = (0.1,0.2,0.0). The 
initial perturbation is Sp — 5.48082 x 10"^ cos(fe-i), Sv — 
2.29279 X 10-6(-2.0, 1.0, 1.0) cos(fc • I), and SB = BoSp. 
The frequency of these shwaves increases as t'^ at late 
times, due to the linear increase with time of both the 
radial wavenumber and va in the presence of a radial 
field. Our algorithm clearly produces convergent results 
on this problem as well. 

^ This corresponds to the maximum growth rate in the magneto- 
rotational instabihty (MRI: IBalbus fc Hawlevlll991h . 



5. SAMPLE NONLINEAR CALCULATION 

The purpose of developing this algorithm is to enable 
new, large shearing box models of disks. Here we describe 
one fruit of this labor: a sample shearing box calculation 
that illustrates the capability of the code. Our model has 
size Lx X Ly X Lz = 8H x 8ttH x 2H.'^ It is unstratified, 
with periodic boundary conditions in the vertical direc- 
tion. The resolution is x Ny x = 128 x 128 x 64. 
The model starts with — (vT5/[327r]) sin(7ra;), so the 
model has zero net vertical field. Velocity perturbations 
of amplitude O.Olcg are added to each zone. 

Figure [T^] shows the evolution of a. As expected, 
the magnetorotational instability grows sharply after a 
few rotation periods. The flow then reaches a nonlin- 
ear regime, followed by saturation. The model saturates 
at a fti 0.01, b roadly consistent with the earlier work of 
iHawlev et al.l ([l995) and others. 

The upper panel of Figure [T51 shows a snapshot of den- 
sity on a two dimensional slice at z = at the end of the 
run t = lOOr^-^. Notice the traihng spiral structures, 
which have an azimuthal extent comparable to the size 
of the box. Also notice that the radial extent of these 
structures is of order H, indicating that, at least in the 
context of this simulation, the correlation function for 
the turbulence is of limited radial extent. We will ex- 
plore this idea further in a later publication. 

For comparison, w e have run the same pro blem with 
the original ZEUS^ (jStone fc Normanl ll992ab). In the 
lower panel of Figure [13] we plot a z = density slice at 
t — lOOfl^^ for the original ZEUS run. Figure [Til shows 
the evolution of the volume averaged magnetic energy 
density for both runs; the run with orbital advection is 
shown as a solid line, while the run with ZEUS is shown 
as a dashed line. There is only a small difference between 
the outcomes visible here, although in the ZEUS run the 
magnetic energy density saturates at a slightly higher 
level. 

One distinct feature of our sample nonlinear calcula- 
tion is the formation of a density dip at center of the box. 
This dip can be clearly seen in the azimuthal and vertical 
averages of the density as a function of x, averaged for 
a period tft — 89.5 — 90.5, as shown in the solid line of 
Figure 1151 In order to improve signal to noise, we time 
average over 11 successive data dumps to generate this 
image. Across the radial grid, the magnitude of the den- 
sity ffuctuation is ~ O.lpo- This density dip has a width 
^ H. Further investigations for the same test problem 
indicate that similar features appear in the results ob- 
tained with other algorithms such as the original ZEUS 
shown in the dashed lin e of Figure [T5|) and ATHENA 
Gardiner fc Stoni [20051 : simulation kindly provided by 
J. Simon). 

This density dip is a generic feature of large shearing 
box calculations. It is associated with large truncation 
errors generated by adverting fluid with respect to the 

^ Most shearing box simulations employ Lx ~ H, although 
some lar ger boxes have been run. Recent examples arc Lx = 
8H dPa paloizou ct al. 2004; Oishi ct al. 2005; Paoaloizou 200i 
IPio ntek fc Ostrikc r ,2007, ). 16-f/ (Kim fc Ostrikcr. ,2006.1 . ITH and 
25H (Kim ct al. 20021 ). although the latter two do not include the 
effects of the MRI. 

^ This version of ZEUS is available at Jim Stone's homepage 
http : //www. astro .princeton. 6du/~ jstone/zeus .html. We use 
this vefsloii in all our code comparfsoh runs. 
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grid. These errors are not distributed evenly in the radial 
direction (Galilean invariance is not satisfied in an Eu- 
lerian integration with shear), therefore large variations 
in density appear in the center of the box because of the 
very small truncation error when w ^ 0. This problem is 
more severe for large boxes without orbital advection be- 
cause the truncation errors increase as x increases. The 
variation in truncation error is relatively small over the 
range |a;| < H/2, so this feature was not observed in 
earlier models with = H . 

For the magnetic field, this also means larger numerical 
diffusivity when |a;| increases. In Figure [TH we plot the 
radial distribution of the spatial averaged magnetic stress 
tensor, time averaged from t — 89.557"^ to t = 90.557^^ 
. 11 data dumps are used to generate this image. Dis- 
sipation of the fields increases with and this leads to 
a gradual decrease of stress (as well as a) towards the 
boundary. At x = ±La;/2 the stress tensor drops to 
~ 50% of its value at the center. 

Notice the strong correlation between the peak of the 
stress tensor and the density dip. This is easy to under- 
stand because in a steady accretion disk aE is a constant, 
where E is the disk surface density. We have also ob- 
served that as the evolution time increases the magnitude 
of the density dip becomes larger due to the accumula- 
tion of truncation errors. For example, from t = AOVL~^ 
to t = lOOrJ^^ the 1st Fourier component ai of the den- 
sity profile ai = J {p/ po) cos{2Trk{x + Lx/2)/ Lx)d'^x in- 
creases from ^ 0.005 to ~ 0.03. In a higher resolution 
study using a, N^x NyX — 256 x 256 x 64 box, oi runs 
from ^ 0.003 to ^ 0.02 over the same time; the feature 
persists, but decreases in magnitude, as the resolution 
increases. 

The radial variation of truncation errors can be seen 
clearly in a linear magnetic field advection test using a 
large, radially extended box. In Figure [17] we plot az- 
imuthally and vertically averaged errors in as a func- 
tion of X for an = lOH box. The alternate appear- 
ances of error minima and maxima are evident. Notice 
that in the orbital advection scheme, numerical errors 
are minimal at those locations x where the relative cell 
shift S ~ —qVLx/S.t//S.y is an integer, because no inter- 
polation is needed. At the box center the fluid does not 
need to be shifted and the errors are minimal; as x in- 
creases, the relative shift gradually increases to 1/2 and 
errors increase to a maximum; beyond this maximum the 
shift then decreases and errors reach a minimum again. 
The ith error minimum should appear at x = Xi which 
satisfies S = —qflxiAt/ Ay = i, where i is an integer. In 
FigurefTTlthe error minima fall exactly at these locations. 

For non-linear large box simulations, one prediction 
for the orbital advection scheme is that the density dip 
should appear at those locations where the cells are 
shifted by an integer amount. In the above = 8H 
model, the relative cell shift S in the orbital advection 
substep is always smaller than one, even at the radial 
boundary. We therefore perform an experiment by ex- 
tending the radial size of the box to 32H. For a size 
Lx X Ly X Lz ^ 32H x 2ttH x 2H box with a resolution 
Nx X Ny X N2 — 1024 X 64 X 64, the estimated time step 
is At ^ 0.0157"^ by assuming \Av\ ^ O.lcs. We then 
estimate that the first integer number shift should occur 
at xi ^ ±7H and the second integer number shift should 



occur at X2 ~ ±14_ff. In Figure [TSl we plot the spatially 
averaged density as a function of x, averaged for a pe- 
riod til = 90 — 100 near the end of the run. The five 
density dips indeed show up at the predicted radial posi- 
tions. It is a coincidence that the locations of the outer 
two density dips are close to the box edge. As shown in 
the above linear advection test, these locations are not 
controlled by the boundary conditions. 

Future large scale shearing box calculations will need 
to eliminate or minimize the numerically induced radial 
variation in mean density. One way of reducing the mag- 
nitude of the dips is to give the whole box a large bulk 
epicyclic motion and let radial oscillations smooth out 
the errors. This may not be an ideal solution because 
any introduced large radial velocity will dramatically de- 
crease the time step. A second approach is to simply 
shift the data by a few H in radius every few f2^^. 

Finally, for larger shearing box simulations our scheme 
is more efficient than the original ZEUS. In the nonlin- 
ear stage of the sample calculation, our scheme is ~ 18 
times faster on a Xeon 3.2GHz machine. Three factors 
contribute to this improved efficiency: (1) By using or- 
bital advection the time step is controlled by Av instead 
of Vorb- The Mach number of the flow with respect to 
a flxed grid at the outer edge of the sample model is 6, 
so the orbital advection scheme reduces the number of 
time steps by ^ 4.8; (2) We implement a larger time 
step than that used in the original ZEUS, which includes 
an unnecessary limit on the time step related to the size 
of the box. This reduces the number of time steps by 
another factor of ~ 2.6. For our sample non-linear calcu- 
lation of size Lj. X Ly X Lz — 8H x 8ttH x 2H box with 
a resolution x Ny x — 256 x 256 x 64, our time 
step is ~ 10 times larger than for ZEUS; (3) We use a 
simpler MOC-CT scheme than ZEUS does, which gives 
an additional factor of 1.3. The remaining factor of 1.1 
is due to minor coding differences. 

6. SUMMARY 

We have developed a scheme for doing orbital advec- 
tion of a magnetized fluid efficiently and accurately using 
interpolation. Our scheme is operator-split, and assumes 
that the magnetic field is discretized on a staggered mesh. 
The main difficulty we have overcome is interpolating the 
magnetic field in a way that preserves V • B = 0. We 
note that other algorithms have been developed for inter- 
polating magnetic fields in a divergence- fr e e man ner. For 
example, iBakard (12Q01D and lToth fc Rod (|2002f ) provide 
prolongation and restriction formulas for interpolating 
fields between grids of different size in Adaptive Mesh 
Refinement codes. Our algorithm is distinct in that it is 
designed specifically for orbital advection. 

Our algorithm can be implemented by encoding the fi- 
nite difference expressions given by equations (|2T|) - ((30| . 
The coefficients in these expressions are defined in Table 
AI. A version of the algorithm implemented in C is avail- 
able at http://rainman.astro.uiuc.edu/codelib.'* 
Our implementation of the algorithm performs orbital 
advection at the end of each time step. In principle, 
however, the orbital advection operator can be inserted 

Note that wo absorb a factor of 1/2 into the definition of the 
van Leer slopes in our code, which introduces a factor of 2 into 
some of the coefficients defined in Table AI. 
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at any point in the series of substeps that make up the 
numerical evolution. We have experimented with differ- 
ent insertion points for some of our tests and have seen 
no significant deviation from the results we present here. 
An important caveat is that the shearing box boundary 
conditions should always be applied at time t before the 
orbital advection substep, and at time t + At after. 

A sample shearing box calculation is shown in Figures 
[TWrSl of the paper. Our scheme produces results en- 
tirely consistent with earlier shearing box calculations, 
but enables the simulation of larger shearing boxes more 
efficiently and more accurately. This should permit the 
study of structures with scales larger than ^ iJ in local 
models of accretion disks. 

One generic feature of these large shearing box sim- 
ulations is the formation of density minima at a; ~ 
in the turbulent stage. We explored the origin of radial 
density variation and have shown that it originates from 
unevenly distributed truncation errors in the radial di- 
rection. Other numbrical algorithms, such as ZEUS and 
ATHENA, are subject to the same numerical artifact. 



The idea b e hind orbital ad vectio ns schemes (see also 
Massed l2000L iGammid 1200 ll and iJohnson fc Gammid 
2005^ is quite general. If (1) the fluid element orbits are 
known at the beginning of the time step (so the interpo- 
lation operator can be constructed), and (2) parts of the 
fluid are moving supersonically with respect to the grid 
(so that orbital advection removes the dominant part of 
the speed that enters the Courant condition) then one 
can in principal obtain a more efficient and more accu- 
rate evolution using orbital advection. 
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of Lawrence Livermore National Security, LLC, (LLNS) 
under Contract No. DE-AC52-07NA27344. This work 
was supported by NSF grant AST 00-03091, NASA grant 
NNG05GO22H, and the David and Lucile Packard Foun- 
dation. C.F.G. thanks the Institute for Advanced Study 
for its support during this work. 



APPENDIX 
A. FLUX CALCULATION 



To flx ideas, consider flrst a given zone of the unsheared grid. For the Cartesian coordinate system we employ here, 
the total flux into the zone is given by 



$ = ^Xijk - ^Xi+ijk + ^Vijk - ^Vij+ik + *%fe - ^Zijk+i, 



where 



ijk 



Ay Az 



1/2 

-1/2 
1/2 



dn„ 



and 



$2;. 



ijk 



-. Ax Az / drix 
J-1/2 

/•1/2 

Aa; Ay / dn^ 

J-1/2 



1/2 

-1/2 
1/2 

-1/2 

1/2 

-1/2 



duz Bx{ny,nz), 



dn^ By{nx,nz), 



duy Bz{nx, Uy) 



(Al) 
(A2) 
(A3) 

(A4) 



The above integrals have been expressed in dimensionless zone units Ux = x/ Ax, Uy = y/ Ay, Uz = zj Az, and the 
integrands are a model for how the field components vary over a zone face. We choose a model that is second-order 
accurate in space: 

Bx{ny, Hz) = bx^jk + bx_dqyijk riy + bx.dqzijk riz, (A5) 

(A6) 



and 



By{nx, Uz) = byijk + by_dqxijk rix + by_dqzijk Uz, 
Bz{nx,ny) = bzijk + bzjdqxi^krix + bz^dqyijkUy, 



(A7) 

where bxijk, by^jk and bzijk are the face-cent ered components of the magnetic field in each zone and, e.g., bx.dqyijk is 
the van Leer slope of bxijk in the y direction (|Stone fc Normanlll992af ). With these definitions for the field components, 
the total flux through a zone face in each orthogonal direction is given by 

^Xijk = bxi-jk Ay Az, ^y^jk = byijk AxAz, ^Zi-jk = bz^jk Ax Ay. (A8) 

Using a subvolume bounded by zone faces from both the sheared grid and the new grid requires, in general, the 
calculation of fluxes through portions of the old grid faces. Figures [2]|4] indicate the subfaces A^;, Ay and A^ over which 
the partial fluxes are defined, and the partial fluxes required for each of the three cases is given below: 



"I>a;TOl = / dydzBx{y,z) = AyAz / duy / duz B 

JAxml Jl/2-nm J~l/2 



^xpl 



dy dz Bx{y, z) — AyAz 



/2 

-1/2+np 
1/2 



dn 



/2 
1/2 



x{ny,nz). 



y I duz Bx{ny,nz), 
-1/2 



(A9) 
(AlO) 



The magnetic field components are defined such that a positive value corresponds to a magnetic fiux into the zone. 
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TABLE 1 

Weight Coefficients 



uiii 



\nm{l-nm) Wig = j^n'^{2np - 3)/s W30 



W12 = Up 



W20 = rim W31 = |nm{l + rim) 

W13 = inp(np - 1) W21 = ^nm(l - rim) 1032 = rip 

-!iii4 = jn^/s ■u;22 = -rip U133 = |np(np - 1) 

i«l5 = gn^C/ - s)/s^ «)23 = -^npCrip + 1) 1034 = -/ 

W16 = j^n^{3 - 2nm) / s W24, = f 11)35 = -^s 

«;i7 = «>25 = -i^s u'36 = ^{f + f + s2/12) 

.1/2 



$zml = / dx dy Bz{x,y) — AxAy / dria; / duy Bz{nx,ny), (All) 

JA,,„i J-1/2 Jl/2-f+n.^s 

I- rl/2 r—l/2-f+n^s 

^zpl= / dx dy Bz{x,y) = AxAy / dna; / driy Bz{nx,ny), (A12) 



/• /-l/Z /•1/2 

$a;m2 = / dydzBx(y,z) = AyAz / dn^ / dnzBx{ny,nz) 



(A13) 



/2-n,„ ^-1/2 
1/2 /.1/2 



$j:p2 = / dy dz Bx{y, z) = AyAz / dnj^ / dn^ Bx{ny,nz), (A14) 
Ja^„2 Ji/2+n„ J-1/2 



/2+np J-1/2 
1/2 (.1/2 



<J>z2 = / dx dy Bz{x,y) = AxAy / dnj; / dnyBz{nx,n. 

Ja,2 J -1/2 Jl/2-f+n^s 



yh 

-1/2 ^l/2-/+n,.s 
-l/2-»im 1-1/2 



(A15) 



$a;m3 = / dydzBx{y,z) = AyAz / dnj^ / driz Bx{ny,nz), (AI6) 

JA^„3 J-1/2 J-1/2 

p |— 1/2+np ^1/2 

$a;p3 = / dy dz Bx{y, z) — AyAz / duj^ / driz Bx{ny,nz), (A17) 

JA^p3 J-1/2 J-1/2 

r rl/2 pl/2-/+n,s 

$28= / dx dy Bz{x,y) — AxAy / dnj^ / duy Bz{nx,ny), (^18) 

Ja^3 J-1/2 J-1/2 

where m and p denote subfaces towards x = —Lx/2 and x = +Lx/2, respectively (to the right and left in Figures [^E]) , 
the fluxes are numbered according to the case for which they are relevant, and the dimensionless zone lengths 

nrn^^ + f , rip EE I - / (Al9) 

are proportional to the azimuthal dimensions of Axm and A^p, respectively.^ 
Using the model defined by equations (jASp through (|A7|) . the partial fluxes are given are given by 

^xmlijk = AyAz {ww bxijk + wn bx_dqyijk) , (A20) 

^xpUjk = AyAz (wi2 bxijk + wi3 bx.dqyijk) , (A21) 

^zmUjk = AxAy (wi4 bzijk + W15 bz.dqx,jk + wie bz.dqyijk) , (A22) 

^zpUjk = AxAy (wi7 bzijk + wis bz.dqxijk + wig bzjiqyijk) , (A23) 

^xm2ijk = AyAz {w2o bxijk + W21 bx_dqyijk) , (A24) 

^xp2ijk = AyAz {w22 bx^jk + W23 bx_dqyijk) , (A25) 

$z2yfe = AxAy {w24 bzijk + W25 bz.dqx^jk + W26 bz.dqyi-jk) , (A26) 

^xmSijk = AyAz (^30 bxijk + W31 bx.dqytjk) , (A27) 

® The origin of the coordinate system for these integrals is defined to be at the location of the field component over which the integral is 
being performed. The integral is over the field in a shea red z one, so that the coordinate axes are parallel to the y, z and sheared x directions 
(the latter axes are indicated by dotted lines in Figures |2I4| I. One can think of an integral over a portion of an x-y subface (e.g., K^ml) in 
the following manner. Imagine the "volume" under the Bz{nx,ny) surface as a series of infinitesimal slabs of length 1 and width dn^ (in 
dimensionless zone units) stacked side-by-side in the radial direction. Integration over riy yields the infinitesimal volume of one of these 
slabs, and a subsequent integration over yields the total volume under the Bz{nx,ny) surface. It is important to perform the integrals 
in the direction of increasing x, y and z so as not to introduce sign errors in the calculation of the fluxes. 



Orbital Advection 9 

^xp3ijk = AyAz {w32 bx^jk + W33 bx_dqyijk) , (A28) 

$z3yfe = AxAy (W34 bzijk + ^35 bz_dqxtjk + W36 bz.dqyiju) , (A29) 

where the coefficients w depend only on the index i (via /) and are defined in Table Al. 

Using Figures as a guide, these definitions can be used to map the sheared grid onto the new grid. The update 
of each magnetic field component can be treated as an independent calculation, although in practice it is natural 
to perform the azimuthal update first, since the updated azimuthal field depends upon the old values for all three 
components. 

Radial Magnetic Field 

The radial flux through a new zone is simply given by the sum of the radial fluxes through the portions of the old 
zones that overlay the new grid. Based upon Figures [2][4l the updated radial flux for each case is given by 
Case 1: 

= '^x^hk - ^xmVljk + '^xmVlj-xk. (A30) 

Case 2: 

= ^<Jk - ^xm2lj, + ^xm2lj^,^, (A31) 

Case 3: 

= <^x^jk - ^xmZl^jt, + $a;m3rj+i,. (A32) 

Converting fluxes to magnetic fleld components via definitions (|A20|) . (|A24p and (IA27P yields the final expressions 
given in the text (equations [2T|-[23]). 

Azimuthal Magnetic Field 

Calculation of the azimuthal field component is the most complicated and requires explicit use of the divergence-free 
constraint. The choice of subvolume over which to sum the fiuxes in a manner consistent with this constraint is not 
unique, so we construct the algorithm under the additional considerations of spatial symmetry and accuracy. 
Case 1: 

We consider three subvolumes for Case 1, indicated by the dark and light shaded regions in Figure [2] and the region 
bounded above and below by Kzmi and A^pi. Summing the fluxes out of the upper (light shaded) subvolume gives 

+ a>zrjfc+i - ^zpl^jk+i - ^-^Jk + '^■^P^uk = 0. (A33) 
Summing the fluxes into the lower (dark shaded) subvolume gives 

+ ^^Ij-u - - ^zpl-^jk+i + ^zplljk = 0. (A34) 

Summing the fluxes out of the region bounded above and below by Azmi and into the region bounded above and below 
by A^pi gives 

+ ^^"^lu-ifc+i - + ^^Plufc - '^^P^ljk+i = 0. (A35) 

Averaging expressions (|A33[) and (|A34p gives the most symmetric algorithm, but the smaller stencil of expression 
(jA35p results in less divergence. The optimum algorithm is therefore to alternate every other time step between the 
average of expressions (R33)) and ((X35|) and the average of expressions (|A34[) and (jA35[) . 
Case 1 (n even): 

= (1/2) ['l^ylj, + ^yZj^,, 2 ^xmllj_,, <i><^, 
+ $a;^ViJfc - 2 '^xplw,, + 2 $zml^V-ifc+i - 2 ^^ml^^^i^ 

+ ^z-j,^, - 2 $ - <i>4jfe + 2 <^zpllj,] . (A36) 

Case 1 (n odd): 

= (1/2) [^yZk + 'i>yTj-ik + ^xIj_,, - 2 ^xmllj_,, - ^x^+,j^,, 

+ ^zZ-ik - 2 <^zmVlj^,, - 2 + 2 ^zpl^j,] . (A37) 

Converting fluxes to magnetic field components via definitions (jA20p - (jA23[) yields the final expressions given in the 
text (equations [24] and [25]). 
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Case 2: 

We consider two subvolumes for Case 2, indicated by the dark and light shaded regions in Figure [H Summing the 
fluxes out of the upper (light shaded) subvolume gives 

+ <i>^2rj_i,+i - $z2rj_u = 0, (A38) 
whereas summing the fluxes into the lower (dark shaded) subvolume gives 

- <i><!7_ife+i + '^z2^j.ik+i + <i><j-ifc - '^zT^j-iu = 0. (A39) 
Taking the average of expressions (|A38p and (|A39p gives 

= (1/2) [^vl,k + <i>2;r/-ifc 

- ^^Ij-ik+i + 2 <i>^2r^_ifc+i) + <^>zlj_^^ - 2 $z2r^_i J . (A40) 

Converting fluxes to magnetic field components via definitions (jA24[) - (IA26[) yields the final expression given in the 
text (equation [IB]). 
Case 3: 

This case, shown in Figure [H is the mirror image of Case 2. Summing the fluxes out of the upper (light shaded) 
subvolume gives 

+ '^^'l+ijk - '^xpi7+iJk + - '^^^Ijk+i - ^^?jk + = 0, (A41) 

whereas summing the fluxes into the lower (dark shaded) subvolume gives 

-$z3i:;fc+i + $z3,'V,. = 0. (A42) 
Taking the average of expressions (jA41[) and (jA42p gives 

^y^r = (1/2) ['fy?jk + <f ya+ifc 

- ^^uk + 2 ^xmS^j, + <i>x^^,j, - 2 <Pxp3t,,j, 
+ ^z^j,^, - 2 $z3l^,fe+i - ^z^j, + 2 $z3rjj . (A43) 

Converting fluxes to magnetic field components via definitions (|A27|) - (IA29[) yields the final expression given in the 
text (equation [17]). 

Vertical Magnetic Field 

The calculation for the vertical field component proceeds in a manner similar to that for the radial component. The 
updated vertical flux for each case is given by 
Case 1: 



$4+1 = $z» , _ <S>zml^j, + <^zml^j_,^ + <i>zpirj+i, - ^zpl^j^, (A44) 



<^zl+' = $4jfe - '^z21j^ + <i>z2,";„i„ (A45) 
<I>z"+i = <^z2jk - $z3rjfe + <^z?,^h+ik ■ (A46) 



Case 2: 
Case 3: 

'^"'ijk 

Converting fluxes to magnetic field components via definitions (|A22p . (|A26p and (|A29p yields the final expressions 
given in the text (equations [2H]-[3D])- 

B. SHEARING WAVE TESTS 

We use the analytical solutions outlined bv lJohnsonI ()2007f ) as the initial conditions for the linear tests in S2](Figures[5]- 
llip . The incompressive solution is given by the real parts of expressions (80)-(82) of that paper. For imaginary uj and 
Co and a Keplerian rotation profile, these are 

6v = 5v cos { k ■ X 



and 



+ 9 (Bl) 
Sva = Sva cos (je ■ X — , (B2) 



with 



and 



where 
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2 2 

6v — jAi [ k , kx ky ^ , kx kz 



k ky 
2ak, 



5vA = --^-^Ai {kl - k^, kxky + 2akl, kxkz - 2akykz) , 



11 
(B3) 
(B4) 

(B5) 
(B6) 

and e is an arbitrary perturbation amplitude. These expressions have been normalized to the correct dimensional 
units. The density perturbation is given by 



Ai = ec.H 



\Lu\n 



n y 2\ijj^\k^ + n^k^z' 



6p Sp ( TT 

— ~ — COS [k ■ X 

Po Po ^ 4 



with 



5p_ 

Pa 



va 5va ^ 20, 

Co Co Cc k 



Sv 



The unstable branch of the incompressive dispersion relation is 

2 



kzfl 



AkvA ■ k 



kM 



and 



[vA ■ ky 



For our choice of initial parameters, va — y^l5/16{fl/kz)z and Hk = 27r(— 2/10, 1/10, 1), these become 



21 



\/67 - 2) ~ lAin^ 



and 



The perturbations in this limit are given by 



16 



1/2 



0.732n. 



Sv = - A(2.)2 



1 



21 1 



21 



/lOl 

VlOO' 50 ' 40a' 5 400a 



and 



with 



and 



Sv, 



16 \uj\ 



VlOO' 50 



2nn \aV67 

1/2 



21 (%/67- 95/16) 
'Vf{VW^^2S~ 



1/2 



0.497. 



(B7) 
(B8) 

(B9) 
(BIO) 
(Bll) 
(B12) 

(B13) 
(B14) 
(B15) 
(B16) 



Dividing through by an overall factor of H^{k^ ~ k^) ~ —(27r)^ (101/100) gives the initial conditions quoted above 
(with e — 10^^ and Cs = O — po = 1). 
We make comparisons based upon the amplitude of the solution, i.e. Sv and Sva rather than Sv and Sva- In 

~ 2 ~ 2 

Figures [51ITU1 then, the quantity that is being plotted is Sv + Sva ■ To extract these quantities from the code, we 
perform spatial sine and cosine Fourier transforms in shearing coordinates on each of the velocity and magnetic field 
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components, and sum the squares of the transforms. As a concrete example, the cosine transform of the radial velocity 
component at time step is given by 

2 

'^^"^^"^ = NNN ^ ^ ^ ''^■'^■'^ ■ ^'^"^ ' ^^^^^ 

^ J* ^ 1=1 j=i k=i 

where k{t) — k{0) + qflkytx. Since the solution as expressed above breaks down as uj transitions from imaginary to 
real, we calculate the analytical amplitudes for the incompressive tests based upon an integration of the full set of 
linear equations/ 

The compressive solution is given by the real part of expressions (83)-(85) of iJohnsonl ()2007l ): 

{Sv, 6va, Sp) = (^Sv, 5vA,Spj cos (k ■ x) , (B18) 

with ^ 

6v ^'^Aci'^k-VA-kvA] , (B19) 



fc2 

, 12 



SvA = ^ A (vA - fc ) , (B20) 



and 

Po 

where 



^ = (iA, (B21) 



y - [va ■ k) cjk^ 

Our choice of initial parameters for this test, va ~ Cs(0.1, 0.2, 0.0) and Hk = 47r(— 2, 1, 1), gives VA-k = Q,so that 
the nonzero solution to the compressive dispersion relation is 

u'^{cl+v\)k\ (B23) 

The perturbations in this limit are given by 

/ ~ \ / I X 1/2 



\^dv,6vA,^^j ^e[vAyA + Pk,VA,l) \^Hk^^^ j , (B24) 
where f3 = c^/va- For our initial conditions {(3 — 20), this is 

which matches the numbers given above (with e = 10^^ and Cg = = po = 1) ■ 
Figure [11] shows the evolution of the azimuthal component of Sva ■ The numerical results are 

Ny , ,n X 

^ k—1 ^ ^ ' ^ 

and the analytical results are calculated within the code using the time dependent version of expression (jB25p . i.e. 

UvAyie]) = VAy^] f Svr./^ \ cos f V uj^'] dtA , (B27) 

V / analytical \ \ ' I \ — / 

\ / \n'=0 / 

with fit" = 0. 

As a final practical consideration, implementing the solutions as described above can introduce divergence into the 
initial conditions. To avoid this, we calculate the vector potential in the Coulomb guage (fc • SA — 0) for the above 
solutions and mmicrically calculate its curl to obtain the initial magnetic field perturbation. The perturbed vector 
potential is 



SA Va ■ k 



A copy of this code is available at http://ralimaii.astro.uiuc.edu/codelib. 



^2a^-l-2a^ + t^)cos(fc-a. + J) (B28) 



Orbital Advection 13 

for the incompressive solution and 

-j== = ^2 sin fc • x) B29 

for the compressive solution. For our initial conditions, these reduce to 



14 \30aVG7 

for the incompressive solution (with a given by expression |B16j ). and 

1/2 



202a, Aa + 105, 40a - y ^ cos (^k ■ x + ^ 



6A c H / 1 \ 

= ( 202a,4a+ 105,40a- ^ ) cos (fc- a; + ^) (B30) 



V4^ 60 \ TT V 14 

for the compressive solution. 



'^^ (2,-l,5)sin(fc.a:) (B31) 
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Fig. 1. — The effect of the baxikground shear flow on a Cartesian grid. The dashed lines represent the old grid after it has been distorted 
by the shear, and the solid lines represent a new grid onto which the sheared grid is to be mapped. 
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Fig. 2. — A slice in the x — y plane for the Case 1 remap. The dashed lines represent the old grid (n) after it has been distorted by 
the shear, and the solid square is a new grid zone (n + 1) onto which the fluxes are to be mapped. The shaded regions correspond to 
subvolumes over which the fluxes are summed for the remap of the azimuthal field. See Appendix |X] for definitions. 
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Fig. 3. — A slice in the x — y plane for the Case 2 remap. The dashed lines represent the old grid (n) after it has been distorted by 
the shear, and the solid square is a new grid zone (n + 1) onto which the fluxes are to be mapped. The shaded regions correspond to 
subvolumes over which the fluxes are summed for the remap of the azimuthal field. See Appendix |X] for definitions. 
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Fig. 4. — A slice in the x — y plane for the Case 3 remap. The dashed lines represent the old grid (n) after it has been distorted by 
the shear, and the solid square is a new grid zone (n + 1) onto which the fluxes are to be mapped. The shaded regions correspond to 
subvolumes over which the fluxes are summed for the remap of the azimuthal field. See Appendix for definitions. 




Fig. 5. — Evolution of the vertical field perturbation for a simple advection test. The thick solid line is the expected result, and the 
thin solid lines correspond to runs at numerical resolutions of 8, 16, 32 and 64 (from bottom to top). The AT = 64 curve is indistinguishable 
from the expected result. 
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Fig. 7. — Convergence as a function of box size with orbital advection on (solid lines) and off (dotted lines). Plotted as a function of 
numerical resolution N is the LI norm of the error in the azimuthal field component with L = H (triangles) and L = lOH (squares). The 
thin solid line is the expected convergence of N~'^. 
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Fig. 8. — Evolution of the radial field perturbation for an incompressive shwave. The thick solid line is the expected result, and the thin 
solid lines correspond to runs at numerical resolutions of Nz = 8, 16, 32 and 64 (from bottom to top). The Nz = 64 curve is indistinguishable 
from the expected result. 
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Fig. 9. — The effects of aliasing for the run shown in Figure |8] for numerical resolutions of = 8 (heavy solid line), 16 (dotted line) 
and 32 (dashed line). The light solid curve is the expected result. 




Fig. 10. — Results from the run shown in Figuro[9]with an additional bulk epicyclic motion superimposed, for numerical resolutions of 
= 8 (heavy solid line), 16 (dotted line) and 32 (dashed line). The light solid curve is the expected result. See text for discussion. 
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Fig. 11. — Evolution of the azimuthal field perturbation for a compressive shwave. The thick solid line is the expected result, and the 
thin solid lines correspond to runs at numerical resolutions of Nz = 8, 16, 32 and 64. The Nz = 64 curve is indistinguishable from the 
expected result. 




Fig. 13. — Density on a z = slice at t = lOOC ^ in the sample nonlinear calculation with orbital advection (upper panel) and ZEUS 
(lower panel). A density dip is visible in both images. 
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Fig. 14. — Evolution of volume-averaged magnetic energy in the sample nonlinear calculation with orbital advection (solid line) and 
ZEUS (dashed Une). 




Fig. 15. — Azimuthal and vertical average of the density as a function of x, averaged from t = 89.50 ^ to t = 90.50 ^ in the sample 
nonlinear calculation with orbital advection (solid line) and ZEUS (dashed line). Both schemes show a density dip at the center of the box. 
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Fig. 16. — Azimuthal and vertical average of the magnetic stress tensor as a function of x, averaged from t = 89.50 ^ to t = 90.50 
in the sample nonlinear calculation with orbital advection. 
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Fig. 17. — Azimuthal and vertical average of the error in Bx as a function of x in a magnetic field advection calculation with orbital 
advection. The error minima appear at those locations where the cell shift in the orbital advection is an integer. 
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Fig. 18. — Azimuthal and vertical average of the density as a function of x near tVl = 90 in an Lx = 32H box with orbital advection. 
The density dips also appear at a; ~ ±7H where the cell shift is ±1, and near the edges where the cell shift is ±2. 



